A seismic simulator used to compute either a forward simulation of a source or the adjoint simulation of a recorded wave field is typically implemented using a time stepping algorithm based upon a selected finite difference approximation to either a first or a second time derivative. Most historical implementations have been for 2nd order time stepping (error is proportional to (Δt)2, where Δt is the time step) because that is easy and efficient to implement and requires fewer resources. Using 2nd order time stepping provides a result with temporal dispersion artifacts. Higher-order approximations are better, because the error for approximation of order n is proportional to (Δt)n which→0 as n→∞ for Δt<1. But any finite-order approximation suffers from some degree of temporal dispersion. If the temporal dispersion is not corrected, the application of this type of simulator for forward simulation or to compute Reverse Time Depth Migration (RTM) images and Full Waveform Inversion (FWI) gradients and Hessians will have errors that degrade the value of these products for petroleum exploration and geophysical prospecting.